Multi-isotopic analysis of zooarchaeological material from Estonia (ca. 200–1800 CE): Variation among food webs and geographical regions

To better comprehend the dietary practices of past populations in the Eastern Baltic region we have created temporally and geographically restricted baselines for the time period of 200–1800 CE. In this multi-isotopic analysis, we report new δ13C, δ15N and δ34S values for 251 faunal bone collagen samples from various archaeological contexts in Estonia representing the most comprehensive set of Iron Age, Medieval and Early Modern Period faunal stable isotope values to date. The results map out the local carbon and nitrogen baselines and define isotopic ranges of local terrestrial, avian and aquatic fauna. We also demonstrate the potential application of sulfur stable isotope analysis in archaeological research. The results demonstrate a clear distinction between δ13C and δ34S values of marine and terrestrial species, however, freshwater fish display notable overlaps with both marine and terrestrial ranges for both δ13C and δ34S values. Herbivores show variation in δ34S values when grouped by region, explained by differences in the local biotopes. This study is the first attempt to connect the Eastern Baltic isotopic baselines and provides more detailed temporal and geographical references to study the local ecologies and interpret the human data.


Introduction
The adoption of stable isotope analysis to study archaeological questions has revolutionized the field by providing new insight into topics as varied as diet, social equality, palaeoecology, and past movements of people and animals. It has become a staple of multidisciplinary investigations, however, some regions such as the eastern part of Europe have been trailing behind in regards to fully embracing the method, as the arrival of modern scientific methods have reached the area with a certain delay. During recent years, there has been a sizeable increase in the number of publications of stable isotope analyses from the Baltic country of Estonia, yet most data have focused on analyzing human skeletal remains [1][2][3][4][5] with little published evidence of related zooarchaeological material [6].

Stable isotope biogeochemistry
The basic principles of stable isotope analysis and its use in archaeology are well established [36,37]. It is based on the premise that bone collagen reflects the isotopic composition of foods (predominantly the protein portion of the diet) consumed over several years before death [38][39][40][41][42]. Carbon isotope ratios (δ 13 C) are used to distinguish between marine and terrestrial sources of carbon, but also between diets based on either C 3 or C 4 plants [36,43]. There is a trophic level enrichment in δ 13 C of about 5‰ between plants and herbivorous animals, and about 1‰ between carnivores and their prey [37,38,44]. Freshwater food webs exhibit widely varying Cisotope signatures, which are generally more negative than marine species but can range in their δ 13 C anywhere between -12‰ and -28‰ [11,45,46].
Nitrogen isotope ratios (δ 15 N) broadly reflect the trophic level of the individual, as δ 15 N values increase by about 3-6‰ with every step up the food chain [36,43,47,48]. Because of the greater number of steps in the marine food chain, fish and marine mammals typically exhibit the highest δ 15 N values. Several factors can influence the nitrogen isotopic composition of plants and animals, e.g., climatic conditions such as precipitation and temperature [49,50], the usage of nitrogen-based fertilizers [51][52][53] and the 'nursing effect' observed in suckling animals [54][55][56][57].
Sulfur isotope ratios (δ 34 S) in terrestrial plants have been shown to vary significantly due to different sulfur sources in the local geology, with additional sulfate sources from the hydrosphere (e.g., ground and stream water) and atmosphere (e.g., sulfur dioxide SO 2 ) all contributing to the average δ 34 S values of the local biologically available sulfur [21,[58][59][60][61]. Freshwater environments can display an equally wide range in their δ 34 S values so that there is a lot of potential overlap between the terrestrial and freshwater ecosystems [18,20,21,59]. In contrast, δ 34 S values of marine organisms living in the ocean are fairly constant and fall within a relatively narrow range between 17‰ and 21‰ [59,62]. Coastal regions can also be affected by the 'sea-spray effect', where the presence of marine sulfates increases soil δ 34 S values, resulting in plants grown on these soils having sulfur isotope values similar to the marine ecosystem [61,63]. This effect decreases with increasing distance from the coast, extending from a few kilometers inland to covering whole islands [63,64].

Overview of sites and samples
A wide selection of sites was included to create a comprehensive isotopic baseline of the whole study area (Fig 1 and Table 1). All distinct regions of Estonia are represented, divided into four sub-regions based on geographical and geological features. 'Northern coastal' and 'western coastal' include sites situated up to 10 km from the (modern) coastline. 'Western coastal' also encompasses the largest islands of Saaremaa, Hiiumaa and Muhu. Inland regions are further divided into 'central inland' and 'southern inland' with the latter roughly corresponding to the geologically distinct area dominated by Devonian sandstone (see Fig 1).
Individual sample selection was based on availability of identified faunal remains in the chosen zooarchaeological assemblages. Altogether 251 samples representing the main sources of animal protein available to the local people were collected, including the most common domestic and wild terrestrial and avian species, and various marine and freshwater fish (Fig 2). Although geese (Anser sp./Branta sp.) specimens are listed under 'wild animals', it is possible that some of these were likely domestic birds [65].
Temporal parameters of each sample are assigned based on archaeological context and, where possible, radiocarbon dating (for 67 out of 251 samples, see S1 Appendix for calibrated dates). To facilitate data analysis, three distinct temporal categories are used (see Table 1), which reflect potential historical differences in settlement pattern, land management, and resource use and availability [66,67]. A best-fit policy was adopted to deal with such a diverse dataset (e.g., when a sample was dated to a transitional period, or the archaeological context spanned a wide temporal range). We recognize that this may obscure the detection of actual temporal patterns, however, this subject can be addressed in future studies. Based on archaeological context, all sites were further categorized as either 'rural' or 'urban'. 'Rural' sites include Medieval and Post-Medieval villages but also all Iron Age sites (e.g., hill forts, settlement and burial sites). The rural societies from the Iron Age until the Early Modern Period were small agrarian villages by their structure that formed larger parishes [67]. In essence, 'rural' sites can be tentatively considered to represent locally raised and caught animals. 'Urban' sites are exclusively from the Medieval and Early Modern periods (including Medieval castles and Hanseatic towns), characterized by higher population density, religious (i.e. Christian) influences, and-in addition to the resources brought from their hinterlandsaccess to more diverse (food) resources, including imported fish.  protocols (with photographs before and after sampling) are stored with the respective collection holders. Specimen numbers and other contextual information can be found in S1 Appendix.

Methodology and quality indicators
The stable isotope and radiocarbon data presented in this paper have been assembled from analytical work conducted in several establishments. Of the total of 251 samples, 197 were processed for collagen extraction in the Biochemistry Laboratory of the School of Natural Sciences and Health at Tallinn University, Estonia. Stable δ 13 C, δ 15 N and δ 34 S isotope analysis for these samples was conducted at the Scottish Universities Environmental Research Centre (SUERC) Radiocarbon Laboratory in East Kilbride, United Kingdom. Radiocarbon ages were obtained from selected samples (n = 13) at the Poznań Radiocarbon Laboratory in Poland. For the remaining 54 samples, collagen extraction and radiocarbon dating was performed at the Leibniz Laboratory for Radiometric Dating and Stable Isotope Research, at the Christian-Albrechts-University of Kiel, and stable isotope measurements done at the Isolab GmbH, Germany. Thirteen (n = 13) samples were further analyzed with Zooarchaeology by Mass Spectrometry (ZooMS) in the BioArCh lab at the University of York, United Kingdom, to differentiate between goat and sheep (S2 Appendix). A full methodological description along with details on the analytical procedure and equipment used is presented in S3 Appendix.
To check sample quality, routinely used criteria were followed [68]. A few samples (n = 8) did not produce enough collagen for δ 13 C, δ 15 N and δ 34 S analysis. Most samples had atomic C: N ratios within the generally acceptable range of 2.9 to 3.6 [68,69], except for two specimens. One of them (cattle-GUsi10837) had a C:N ratio of 3.7 but satisfactory elemental concentrations, so was included in the dataset. The other (hare-GUsi10812) had a ratio of 3.8 but since this sample also had problematic δ 34 S quality indicators, it was discarded from further analysis. For carbon and nitrogen concentrations, most samples were above 30% for %C and more than 10% for %N, indicative of good quality collagen [68]. Seven samples had both values below that range but still within the accepted lower limits of 13% for %C and 5% for %N [70].
Quality indicators for sulfur stable isotopes are not as established as for δ 13 C and δ 15 N. Arguably, the research concerning these criteria is still in progress as the number of controlled, large-scale studies on δ 34 S isotopes is relatively low. In our study, mammalian and avian samples have an average C:S atomic ratio of 482 ± 96 (1-sigma standard deviation), N:S atomic ratio of 144 ± 29 (1SD), and sulfur concentration (%S) of 0.24% ± 0.1% (1SD), consistent with results reported by other researchers [14,17,23,71]. We have excluded nine samples with C:S/ N:S ratios of <250/75 and %S >0.40% as these values may represent sulfur contamination. The δ 13 C and δ 15 N values for these samples were not excluded as the other collagen quality indicators were acceptable.
Fish atomic C:S and N:S ratios reported here (on average 254 ± 34 and 81 ± 12, respectively) are somewhat higher than those seen in other studies [18,23,71] and the %S lower (0.42% ± 0.1%). However, all fish samples have acceptable %C and %N content, C:N ratios and visual appearance. There has not been enough research done on fish δ 34 S to argue whether the reported values may be influenced by contamination. Since it seemed to us that these data were still meaningful, we chose to include them in our study.
Statistical analysis was applied to determine whether differences between groups of data were (statistically) significant. A non-parametric Mann-Whitney U test was preferred due to the presence of outliers and non-normally distributed data (as assessed by the Shapiro-Wilk test of normality). When more than two groups were compared, the non-parametric Kruskal-Wallis H test was applied. The rejection of the null hypothesis for the Kruskal-Wallis test was automatically followed by a post hoc analysis using Dunn's procedure with a Bonferroni correction for multiple comparisons. Statistical uncertainty was set at 95%.

Results and discussion
Results of stable isotope analyses are summarized in Table 2 and plotted in Figs 3-5. A full list of analytical measurements along with associated radiocarbon ages, contextual and collagen quality information can be found in S1 Appendix. Additional summary statistics and results of statistical tests of significance are presented in S4 Appendix.
Due to the diverse dataset, it was not possible to follow a rigid sampling procedure, e.g., only sampling specific skeletal element(s) or only adult (mature) animals. We acknowledge that these factors may result in data that are not entirely comparable due to minor isotopic differences caused by intraskeletal and interspecies variation in bone collagen turnover rates [72][73][74][75] and age-specific diets [54][55][56][57]. However, the discrepancies in laboratory protocol and analytical equipment between the two sets of analyzed data should not hinder their meaningful comparison as has been demonstrated by various inter-laboratory experiments [76,77]. It is unlikely that the listed differences are significant enough to obscure any major patterns that may exist in the local baseline. Table 2. Summary statistics of stable carbon, nitrogen and sulfur isotope data of analyzed faunal samples by group (number of samples that produced results, average, 1-sigma standard deviation, minimum and maximum values, range). Goose δ 13 C and δ 15 N results have been previously published [60].

PLOS ONE
5.1 ± 1.7‰, and geese -22.7 ± 0.3‰ and 9.9 ± 2.1‰ (Table 2 and Fig 3). These results reflect a temperate, C 3 -plant dominated ecosystem, characteristic to local flora. Domestic herbivores and omnivores have statistically significantly different δ 13 C and δ 15 N values (Mann-Whitney U test, p<0.03 for both variables) due to the inclusion of animal protein in the diet of omnivorous animals. Domestic herbivores. Among domestic herbivores, there is a particularly large variation in both δ 13 C and δ 15 N values. δ 13 C values range from -20.1‰ to -23.1‰. There are no statistical differences in δ 13 C values between domestic herbivores found from urban vs. rural contexts (S4 Appendix). Lower δ 13 C values (below -22‰) were likely due to forested environments [78,79], e.g., animals left to graze freely in the surrounding woodlands, whereas higher δ 13 C values can be associated with animals raised in open areas, probably close to human settlements. However, there are no statistically significant differences between domestic herbivore δ 13 C values by region, nor when only comparing coastal vs. inland samples. Although terrestrial plants display geographical variation in stable carbon isotope ratios due to climatic conditions such as water availability, temperature, and sunlight exposure [40,80,81], these factors probably influence herbivore isotope values on a much larger scale than observed here. Herbivore carbon isotopic values also show no statistical differences between periods, suggesting relative consistency in herding conditions from the later Iron Age to Medieval and Post-Medieval times.
Domestic herbivore δ 15 N values range from 2.7‰ to 10.3‰. The large variation can be explained by various scenarios which are not necessarily mutually exclusive. For example, atypically high δ 15 N values in herbivores can be caused by nursing, i.e. suckling animals literally feeding off their mothers [54][55][56][57]. Regrettably, it is not always possible to determine the exact age of the animal (i.e. whether it was a suckling) based on a single bone or bone fragment. Of the two herbivores that can be considered as outliers based on their high δ 15 N (and The usage of traditional fertilizers such as animal manure on plants can also increase their δ 15 N values by a full trophic level [51][52][53]. If this enrichment is carried up the food chain, it can lead to herbivores having δ 15 N values characteristic to omni-or carnivores. Low δ 15 N values could thus reflect consumption of wild plants, whereas higher values can be due to animals grazing on cultivated land or being directly fed manured crops. Statistical analyses demonstrated that urban domestic herbivore δ 15 N values were significantly different from rural herbivore δ 15 N values (Mann-Whitney U test, U = 1468.5, p = 0.007) with the latter having on average higher δ 15 N values (Table 3). This would imply that rural herbivores were consuming more intensively manured plants or that they were more often slaughtered at young age compared to urban herbivores. However, the distinction between urban and rural domestic herbivores may be altogether misleading as even in the largest towns, herbivores usually grazed outside the town limits, and livestock may have been brought to towns for butchering or selling at the market from the surrounding rural areas.
In A general trend in global soil δ 15 N values has been recognized, which seems to decrease with increasing annual precipitation and decreasing mean annual temperature [49]. Similarly, European terrestrial herbivore δ 15 N values show a distinct decline during the Last Glacial Maximum compared to earlier and later phases when the climate was milder [50]. The only significant climatic fluctuation during the study period was the Little Ice Age [82,83], which could have theoretically resulted in decreased herbivore δ 15 N values. However, the LIA peaked in the Early Modern period (ca. 16-19 th centuries), having been preceded by relatively warmer conditions in Northern Europe during the so-called Medieval Warm Period (ca. 9-14 th centuries)  [84]. As such, we would expect this to have produced an opposite δ 15 N trend to what is observed here. It thus seems unlikely that climate had a significant effect on terrestrial herbivore nitrogen isotopic values in our dataset.
Considering that many of the included sites were in continuous use, non-radiocarbon dated samples may have been affected by potential mixing of earlier and later material. However, only comparing radiocarbon dated samples does not significantly change the outcomes regarding average group values and results of statistical tests of significance, nor are the results notably affected when removing the two outliers (S4 Appendix). There are also no consistent differences between Iron Age and Medieval δ 15 N values when analyzed by species (see S2 Fig,  S4 Appendix). This taken together implies that there are no systematic differences in domestic herbivore nitrogen isotope values between Iron Age and Medieval periods, and that the underrepresentation of Iron Age (and Early Modern) samples may have skewed the results of statistical comparisons.
Some of the variation in domestic herbivore δ 13 C and δ 15 N values is certainly due to interspecies differences, although statistically only Medieval period goats and sheep had significantly different δ 15 N values (S4 Appendix; all other comparisons between species, either combined or by period, failed to produce significant differences). Additionally, mean values of herbivorous species are not that different from each other, and it is evident from Fig 4 that intra-species isotopic variation is sometimes even greater than inter-species variation, suggesting that species-specific differences were not the main factor behind the observed variation. The wide variation in domestic herbivore δ 13 C and δ 15 N values may have also been due to sampling different sites across a large area. The notable range of herbivore isotopic variation seems to indicate that domesticates were handled in varying conditions depending on available resources and conditions specific to any given site.
Domestic omnivores. For domestic omnivores, all three species-the pig, dog and chicken-are distinct in their isotopic distribution pattern. The wide range of isotopic variation for omnivores is to be expected considering they have the option to alternate between herbivorous and carnivorous diets. δ 13 C values range from -17.9‰ to -22.8‰ (4.9‰ total) and δ 15 N values from 5.3‰ to 11.7‰ (6.4‰ total). The range of omnivore δ 15 N values is in fact smaller than was reported for domestic herbivores (7.6‰; although this was due to the inclusion of two outliers with high δ 15 N). There are no statistical differences between urban and rural omnivore δ 13 C and δ 15 N values. Statistical comparisons between omnivores from different regions and periods were run (S4 Appendix), but due to the limited number of samples in some of the groups, the results may be skewed and will thus not be discussed further here.
Pigs show the most variation in δ 15 N values, and dogs in δ 13 C values. Pigs also have the lowest mean δ 13 C and δ 15 N values and dogs the highest, whereas chickens (although with a limited number of samples) tend to plot in-between. There is a positive significant correlation between pig δ 13 C and δ 15 N values (Pearson's R = 0.43, p = 0.011 for all values, but R = 0.77, p<0.001 when removing four outliers: GUsi9209, KIA-55257, GUsi10831, GUsi10844), suggesting that variation in pig isotope values can be best explained by increased reliance on higher trophic level protein. Pigs sampled here range from almost fully herbivorous (possibly left to roam freely in forests) to those similar to chicken and dog. However, most pigs show little overlap with ruminants when compared with samples from the same site (e.g., urban areas such as Tallinn, Tartu, Viljandi and Pärnu where at least ten domesticated animals were analyzed), suggesting that those on predominantly herbivorous diets were not consuming the same plants as (most) cattle and caprids were.
Dogs are the most distinct of all omnivores, indicating a unique diet, and are the only terrestrial species with δ 13 C values above -20‰. Dogs' isotopic values can sometimes be used as a proxy for human diet due to their role as companions and pets. Indeed, published human isotopic data from two historic period sites in northern Estonia [5] exhibit notable resemblance to dogs' values, especially for urban contexts. Individuals buried at the suburban St Barbara cemetery (15 th -18 th century) have mean values of -20.0‰ and 11.2‰ for δ 13 C and δ 15 N compared to -20.0‰ and 10.7‰ for urban dogs in this study. For rural dogs, the mean values of -20.5‰ and 9.6‰ for δ 13 C and δ 15 N are more distinct from rural humans buried at the Kaberla village cemetery (12 th -17 th century), who have mean values of -19.8‰ and 10.4‰. The larger discrepancy between rural dogs and humans may be due to the fact that the cemetery is from northern coastal Estonia whereas the majority of rural dogs sampled for this study are from inland sites further south. However, the close similarity of urban commoners' and dogs' isotope values from similar contexts supports the assumption that dogs were likely fed household scraps and leftovers, including some aquatic resources as evidenced by two dogs (both from the same urban coastal Medieval context) with δ 13 C values of -17.9‰ and -18.3‰.
Wild mammals. Wild mammals have, in general, lower δ 13 C and δ 15 N values compared to most domesticated animals. Statistical comparisons between wild and domestic herbivores demonstrate significant differences in their δ 13 C values (Mann-Whitney U test, U = 149, p<0.0001), but not for nitrogen (U = 383, p = 0.066). No similar comparisons were run for omnivores due to their fewer numbers in the dataset. The isotopic range of wild mammals reported here is characteristic of forested environments and the consumption of wild (i.e. non-cultivated) plants. There is a noticeable range of isotopic variation observed for hares, most of which are from similar contexts (Medieval southern inland), probably reflecting differences in habitat and food availability. One hare is a clear outlier (GUsi10234, see also Fig 4) that may have been kept in captivity (e.g., as a pet).
Wild boars have isotope values that partly overlap with some of the domestic pigs, which supports the assumption that pigs were sometimes left to roam freely in forests. While it is often not possible to distinguish between the wild boar and domestic pig in the zooarchaeological material, the animals from the periods under study are quite distinctive in their size, with Iron Age and Medieval pigs being much smaller than their wild ancestors.
Geese. Geese comprise a separate group as domestic and wild geese are hard to differentiate based on bone morphometrics alone. While at least some of the geese sampled here were most likely domestic [65], specifically those with high δ 15 N values indicative of fattening, it is still noteworthy that geese occupy a completely unique niche among other sampled terrestrial fauna. Their low δ 13 C values are indicative of influences from freshwater environments and/or uncultivated landscapes.

Sulfur isotopic variation in the terrestrial food web
Terrestrial animals display a remarkably wide range (17.8‰) in their δ 34 S values, from -6.3‰ to 11.5‰ (average 4.8 ± 3.8‰) (Figs 5 and S3). There is also a great deal of overlap between different species. Terrestrial animals as a group have statistically significantly different δ 34 S values compared to both marine and freshwater species (p<0.0001).
Animals obtain their sulfur through dietary protein (either plant-or animal-based) with only a minor (on average +0.5‰) trophic level change between consumer and food source [44,61,85]. Therefore, animals eating in a specific region should have δ 34 S values similar to the surrounding terrestrial vegetation. Because of this and how the local sulfur isotopic baseline has remained relatively constant in (pre)history, the observed variation in our dataset is unlikely to have been caused by species-or period-specific differences. For this reason, sulfur isotopic data should be analyzed by region.
In Table 4, average terrestrial δ 34 S values from the four regions are compared with each other, demonstrating noticeable and statistically significant (except for between 'northern coastal'-'western coastal' and 'western coastal'-'central inland') differences between them (see S4 Appendix for results of statistical tests). There appears to be a north-west/south-east gradient with δ 34 S values increasing further south, ranging from -6.3‰ in 'northern coastal' to 11.5‰ in 'southern inland'.
Since the terrestrial fauna also includes omnivorous animals, whose δ 34 S values may have been influenced by direct or indirect consumption of aquatic protein, the regional average values may not display the true terrestrial sulfur signal obtained through local plants. In Table 4, δ 34 S values of all terrestrial species are further compared with only herbivorous animals (both wild and domestic) and only omnivorous animals (including geese). While geese are predominantly herbivorous, their δ 34 S values cannot be considered a reliable indicator of the local terrestrial baseline due to their preference for aquatic plant consumption, in addition to the possibility that some sampled geese may have been wild (i.e. migratory) birds. However, removing the geese samples from analysis altogether does not significantly change the outcome, as most specimens actually plot near the regional mean δ 34 S value, which may suggest that they were locally raised.
The differences between regional average δ 34 S values remain consistent even if only terrestrial herbivores are included in the analysis (Table 4). While the mean values and 1SD are very similar between the 'all terrestrial' and 'herbivore' groups, again, statistically 'northern coastal'-'western coastal' and 'western coastal'-'central inland' are not significantly different.
Omnivore mean δ 34 S values display a slightly different pattern with coastal regions demonstrating distinct (but not statistically significant) differences between omnivore and herbivore mean δ 34 S values, with the former having on average around 2‰ higher values compared to the latter. This variance is best explained by the consumption of aquatic protein (with higher δ 34 S values) by omnivorous animals, which has resulted in an increase of local terrestrial δ 34 S values. Inland regions do not show these differences. It is unlikely that this is due to inland omnivores not consuming any aquatic protein; the lack of difference is rather caused by the local terrestrial δ 34 S baseline already being higher there.
The increase in average coastal omnivore δ 34 S values may also reflect an inclusion of dietary sources affected by the 'sea-spray effect', which results in increased δ 34 S values due to the effects of airborne marine sulfates on coastal soils and plants [61,63]. However, no consistent sea-spray effect was observed among coastal herbivores. Additionally, coastal fauna have on average the lowest δ 34 S values in the whole dataset, although some herbivores from coastal areas do have noticeably high δ 34 S values. The absence of a consistent sea-spray effect can be taken to suggest that either its range is relatively limited and/or that most of the sampled terrestrial fauna were not feeding close to the shoreline. However, without sampling local geology, the true extent of the sea-spray effect is difficult to estimate. The terrestrial baseline may theoretically be highly negative but dampened by various proportions of marine sulfates.
Looking at individual terrestrial faunal δ 34 S values by region (Fig 6), 'northern coastal' displays the largest range of values (16.5‰), from -6.3‰ to 10.2‰. Most herbivores have δ 34 S values within the 1SD range of the regional average, however, many samples fall outside this range, including an elk (with a δ 34 S value of 6.6‰) and both herbivorous and omnivorous domesticates. It is noteworthy that two dogs with clear signals of marine resource consumption (with δ 15 N values around 12‰ and δ 13 C values around -18‰) also have some of the highest δ 34 S values in this region. This can be compared to another dog sample from the same context that displays much more 'terrestrial' isotopic values (9.2‰ and -20.6‰ for δ 15 N and δ 13 C, respectively), coupled with a very low δ 34 S value of -2.2‰. Based on the size of the sampled bone, this individual may have also been a wolf. It would thus seem that the lower end of δ 34 S values reported for 'northern coastal' may actually represent the local terrestrial baseline, whereas higher δ 34 S values can be explained by the influences of the sea-spray effect (for herbivores and omnivores) and consumption of aquatic resources (for omnivores). Terrestrial fauna from the 'western coastal' region displays equally marked variation in δ 34 S values, albeit with a slightly smaller range compared to the northern coast -13.8‰, from -2.9‰ to 10.9‰. There are also notably less samples that fall outside the 1SD range for 'western coastal'. A few domestic herbivores have δ 34 S values below 0‰. These animals may have been herded in areas closer to the northern coastline. On the other hand, herbivores with the highest δ 34 S values were likely raised on or near coastal soils influenced by the sea-spray effect (although it is also possible these animals originated from more southern/inland parts of the land). Interestingly, several geese (sampled from a single context) all have very similar δ 34 S values (furthermore, their δ 13 C and δ 15 N values are also very much alike). These particular geese were probably domestic and kept for food since their high δ 15 N values and large bone size were indicative of fattening [65]. It seems likely that these birds were kept locally and possibly on a controlled diet (i.e. not left to graze freely) so that their uniform δ 34 S values may be considered a good representation of the local baseline (possibly with a small enrichment if these geese were also consuming aquatic plants or fed fish offal).
In 'central inland', δ 34 S values again show a rather substantial range (11.4‰, from -1.5‰ to 9.9‰) but the results are much more concentrated compared to coastal regions. A few animals with values below 0‰ are likely from areas closer to the coast, and likewise samples with very high δ 34 S values may originate from other areas (either from the coast where they were affected by sea-spray, or from further south).
For the 'southern inland' region, the range of δ 34 S values is narrow compared to other areas (8.6‰, from 2.9‰ to 11.5‰). Several wild hares, which likely were local specimens, were sampled here and most fall within the 1SD range of the regional average, including the outlier with atypical carbon and nitrogen isotope values. A few herbivores have δ 34 S values around or above 10‰ which may suggest that the local terrestrial baseline might be even more positive further south/southeast. The highest reported δ 34 S value (of any terrestrial animal in this study) is from an elk (11.5‰). On the other hand, a beaver has a rather low δ 34 S value of 3.7‰, similar to two geese and a few other samples with δ 34 S values below 5‰. Fauna with such low δ 34 S values was probably non-local, brought to the 'southern inland' region through trade or hunting activities.
In baseline of a certain region, it is clear from our data that geology seems to be a major factor behind the observed sulfur isotopic variation in this region. Further studies on rock, soil and plant δ 34 S values would be beneficial in supporting this hypothesis.
The fact that there exists such a gradient in local terrestrial faunal sulfur isotope ratios from the north-west to south-east can be extremely useful for studies of human migration and

PLOS ONE
animal trade. However, this assumes that we understand all the factors influencing terrestrial faunal sulfur isotope values and can determine the local baseline with as much accuracy as possible.
The exact range of δ 34 S terrestrial baselines in coastal regions requires additional research, as the current variation among sampled herbivores is far too great and cannot be reliably used as an indicator for sulfur isotope values of local plants. However, it is reasonable to assume that the lowest reported δ 34 S values from coastal regions likely reflect the local terrestrial baseline and the remaining variation can be explained by marine influences (i.e. sea-spray). We do not know the true extent of the sea-spray effect, but in the 'northern coastal' region most animals with δ 34 S values below 0‰ originate from sites to the south and east of Tallinn which can be considered to represent the rural hinterland, whereas herbivores from the urban center of Tallinn have an average value of 3.3‰. In fact, all samples from the 'northern coastal' region with δ 34 S values above 6‰ are from Tallinn, a town situated right on the coastline.
A similar situation can be observed in the 'western coastal' region where most fauna with high (above 6‰) δ 34 S values are from another coastal town -Pärnu. It is very likely that livestock from Tallinn and Pärnu were sometimes herded on coastal meadows (a practice that is also known from recent historical periods [86,87]) which would have been enriched with marine sulfates.
In the context of using sulfur isotope analysis to reconstruct diet and origin of humans from these periods, average δ 34 S values specific to towns and smaller regions shown in Fig 8 can be used as a reference. The local animal range (i.e. the terrestrial baseline) can be defined based on their average δ 34 S values ± 1SD. The 1-sigma standard deviation is much higher for coastal fauna but we do not believe this needs correcting since animals herded on coastal soils (with enriched δ 34 S values) were also routinely consumed and thus the larger variation would also be reflected in local human δ 34 S values. This has significant implications for using δ 34 S for palaeodietary reconstructions in this region, especially in terms of marine resource consumption, as it would be difficult to fully distinguish whether increased δ 34 S values of humans are due to consumption of aquatic resources or terrestrial resources influenced by the sea-spray effect.
For samples from 'central inland' region, a distinction should be made between northern and southern sites (broadly corresponding with the Ordovician/Silurian boundary shown in Fig 8) when using faunal data for local baseline reference. Northern central sites have lower δ 34 S values (on average 3.7‰) compared to southern central sites (6.0‰), although the difference is not statistically significant (Mann-Whitney U test, U = 70.5, p = 0.082).
The fact that terrestrial fauna in the 'southern inland' region has high δ 34 S values approaching those seen in freshwater and marine environments, means that this proxy cannot be used to detect aquatic resource consumption in the region. However, non-local animals and migrants from northern and western regions are much more likely to be identified here, as their low δ 34 S values will appear as outliers.

Isotopic variation of aquatic species in and around the Baltic Sea
The results of this study display significant variation in δ 13 C, δ 15 N and δ 34 S values of aquatic species, with seals having average values of -16.1 ± 0.5‰, 13.6 ± 1.1‰ and 16.5 ± 1.2‰, marine fish -14.5 ± 1.4‰, 12.5 ± 2.4‰ and 15.9 ± 1.5‰, and freshwater fish -18.6 ± 5.2‰, 10.4 ± 2.0‰, and 12.4 ± 2.0‰, respectively (see Table 2 and Figs 3 and 5). Marine and freshwater species are statistically significantly different in their mean δ 15 N and δ 34 S values (p<0.004 for both), but not for δ 13 C (S4 Appendix). Such a wide variation suggests that fish from a wide range of habitats were exploited for food. These potential habitats are indicated in Fig 9, where δ 13 C and δ 15 N values of analyzed species group to form distinct niches.
Five cod samples (all of them from large sized individuals) are clearly separate from the rest of the fish, displaying a 'typical' marine signal with δ 13 C values around -14‰ and δ 15 N values around 14‰. The combination of their size and isotope values indicates that these fish were most likely imported from the North Sea [16]. The remaining cod samples have highly variable δ 13 C values (from -14‰ to -17‰) coupled with relatively low δ 15 N values (ca. 10‰). We presume these fish to be caught locally, i.e. to represent the isotopic signal of the Baltic Sea. This would be in accordance with previously published data [16,[88][89][90] that showed Medieval cod originating from the eastern Baltic (specifically from sites in Poland and eastern Sweden) have δ 13 C values around -16‰ to -18‰, whereas cod from western Baltic and the North Sea have carbon isotopic values around -13‰ to -15‰. Since carbon isotope ratios are directly correlated with water salinity [32], various marine consumers from the brackish Baltic Sea have been previously shown to have more negative δ 13 C values compared to 'typical' marine environments [14,15,91].
Cod caught from the Baltic Sea vs. from the North Sea are distinguished by both their δ 13 C and δ 15 N values. While the former variation is due to the brackish conditions (i.e. decreased salinity) of the Baltic Sea, lower δ 15 N values of locally caught cod may be an indication of their trophic position. North Sea cod have higher δ 15 N values because of the increased number of trophic levels in the oceanic food web. They also have faster growth rate, resulting in oceanic cod being generally larger in size [92,93]. Since the size of the fish is also closely associated with trophic levels, it is thus correlated with δ 15 N values [94]. Among the four cod specimens from the Baltic Sea, two are identified as medium-sized and two as large (when compared with local reference individuals), even though they all have relatively uniform δ 15 N values.
Two marine samples in our study (a cod and a turbot) display a combination of low δ 15 N values (resembling cod from the Baltic Sea) and high δ 13 C values (similar to cod from the North Sea). In these instances, the δ 13 C values probably do not reflect variations in seawater salinity, but rather changes in water turbidity. For example, δ 13 C of aquatic species varies based on whether they are from faster moving open water areas (offshore) or from slower moving waters (nearshore), so that very high δ 13 C values are often seen in fish living in stale waters [94]. Turbot, in fact, is a fish that prefers shallow waters near the shore [95], whereas cod are also known to sometimes favor coastal, benthic habitats [96].
Seals sampled in our study have generally uniform isotope values, suggestive of feeding on local (i.e. eastern Baltic) marine fish, except for the harp seal from Rebala (KIA-55657) which has higher δ 13 C and δ 15 N values (-15.4‰ and 14.9‰, respectively) compared to the other three seals (gray and ringed seals). This range of values is similar to those of imported cod and is somewhat higher than reported for Bronze and Iron Age harp seals that were feeding and breeding in the Baltic Sea [97]. However, the differences are not that significant to claim that our harp seal was necessarily non-local (i.e. migratory).
For freshwater fish, which are characterized by even greater variation in both their δ 13 C and δ 15 N values as compared to marine species, two distinct niches are clearly visible in Fig 8. Carbon isotopic variation in freshwater fish is common for species that are flexible in their occupying niches and that can live in both fast flowing riverine, nearshore, and offshore waters [94]. Variability in fish δ 15 N values, however, can be influenced by fish size but also by environmental factors such as oxygenation, pH, and nitrogen limitation (e.g., influences of wetland and upstream environments, human activities, etc.) [94].
One of the two niches (including pike, perch and bream) is characterized by carbon isotopic values typically associated with freshwater fish, with low δ 13 C values (from -23‰ to -25‰) and δ 15 N values between 10‰ and 14‰. These fish probably originated from fast flowing rivers in inland regions. The noticeable variation in δ 15 N values even within one species (e.g., pike or perch) may reflect individual differences in size. While fish samples that were identified as having belonged to large individuals had higher than average δ 15 N values, high nitrogen isotope values were also seen in those specimens that were identified as medium-sized. Most freshwater species with low carbon isotope values are from inland sites.
The other group of freshwater fish is set apart by their high δ 13 C values (up to -11‰). These include pike, perch, bream and ide. The carbon and nitrogen isotope values for this group partially overlap with those from locally caught cod (the two groups are statistically similar, p>0.08 for both), suggesting that these samples represent freshwater fish living in the brackish conditions of the Baltic Sea. The highest δ 13 C values are probably influenced by a littoral habitat, i.e. living in shallow waters close to the shore, as also suggested by Eriksson et al. [15], who reported very high δ 13 C (between -13‰ and -10‰) values for pike and perch from the Baltic Sea island of Ö land (Sweden). The low δ 15 N values of species such as bream and ide reflect their feeding habits and lower trophic position compared to piscivorous species such as pike and perch. Most freshwater species with carbon isotope values higher than -20‰ are from coastal sites. Figs 10 and 11 show δ 34 S values of aquatic species plotted against their respective δ 13 C and δ 15 N values. Marine species in our study are well distinguished by their δ 34 S values, generally falling between 15‰ and 18‰. This range of δ 34 S values was also measured for modern cod caught from the Baltic Sea [98]. A cutoff point of 15‰ is generally used for δ 34 S to distinguish between fully marine and coastal environments, with another cutoff point of 10‰ to signify marine environments with substantial riverine freshwater influx such as estuaries [14,[98][99][100].
Interestingly, variation in δ 13 C and δ 15 N values between cod from the Baltic Sea and North Sea is not reflected in their δ 34 S values. While this may be an effect of the small sample size, it seems that marine δ 34 S values are not correlated with salinity the same way that δ 13 C values are, or that there are other factors behind the observed variability. The latter seems probable as even among cod thought to be imported from the North Sea or beyond (with relatively

PLOS ONE
uniform δ 13 C and δ 15 N values) there is noticeable variation in δ 34 S values, from 12.2‰ to 17.6‰. The Rebala harp seal with isotope values similar to imported cod also has a δ 34 S value not reflective of the 'typical' marine sulfur isotope range (14.8‰), while the other seals believed to be local have much higher values (from 16.8‰ to 17.4‰).
The apparent overlap in oceanic vs. Baltic Sea δ 34 S values is all the more surprising considering that local freshwater fish living in brackish conditions have statistically significantly different δ 34 S values (but statistically similar δ 13 C) compared to purely marine species from the Baltic. While this may again be affected by the size of our dataset, the majority of brackish freshwater fish have δ 34 S values in the range of 10‰ to 15‰. This suggests that the freshwater fish mostly occupied coastal waters whereas cod and seals preferred offshore areas unaffected by freshwater sulfates. This may also explain why some of the specimens most depleted in 34 S are also the most enriched in δ 13 C, as both of these instances occur in near-shore waters which likely experienced significant mixing of marine sulfates with freshwater and terrestrial sulfur sources.
When Leakey et al. [100] compared δ 34 S values of fish caught from the Thames estuary to those from adjacent coastal regions in southeast England, they observed a clear gradient of increasing δ 13 C and δ 34 S values that resulted from the mixing of marine and freshwater, with δ 34 S values ranging from -5‰ to 15‰ across the gradient. In our dataset, no similar gradient is evident when comparing freshwater species from purely freshwater (e.g., riverine) to brackish (coastal) waters. In fact, 'freshwater' and 'brackish' fish show significant overlap and no statistical differences in their respective δ 34 S values (Figs 10 and 11).
On the other hand, freshwater fish from inland locations have been shown to have highly variable δ 34 S values from site-to-site, even within one region [18,20,21]. It is possible that the observed similarity between purely freshwater and brackish fish δ 34 S values is coincidental and that the samples reported here reflect local riverine sulfates of various bodies of water. More research on fish from well-defined, single contexts in Estonia is necessary to further explore the variation in the local aquatic baselines, as there is a distinct lack of paired δ 13 C, δ 15 N and δ 34 S isotope data from the Baltic region.

Conclusions
The large dataset analyzed here has established faunal baseline ranges to better understand the reasons behind the observed large variation in isotope values and to aid in future (human) paleodietary research of material from related contexts. Our results demonstrate that domestic omnivores and herbivores are well distinguished by their δ 13 C and δ 15 N values, although there is a great deal of intraspecies variability. No regional trends in carbon and nitrogen isotope values are evident, however, the significant variation in domestic herbivore δ 15 N values require further research to determine the exact extent and significance of the observed pattern.
This study has demonstrated the importance of sulfur stable isotope analysis in local archaeological research by showing clear differences in marine vs. terrestrial δ 34 S values and the existence of a north-west/south-east gradient in terrestrial δ 34 S values influenced by local geology. Low terrestrial δ 34 S values (ca. 0‰) can be interpreted as originating from coastal regions (especially the northern part of the land), whereas highly positive δ 34 S values can be caused by various factors such as aquatic resource consumption, marine sulfate influence on coastal soils, or originating from the southern inland regions. These differences can potentially be exploited to study migration at a local level, e.g., to shed more light on exogamous practices that were common in (pre)history.
Although some researchers have successfully used sulfur isotope analyses to identify freshwater resource consumption, this is problematic in our study area due to δ 34 S values of freshwater fish partly overlapping with both terrestrial and marine sulfur isotopic ranges. However, δ 13 C values of freshwater species seem to be a good indicator of their habitat, e.g., whether they originated from inland (e.g., riverine) or coastal (e.g., estuary) regions. Since aquatic resources with such diverse isotopic signals were all routinely consumed, identifying their respective signals in consumer (i.e. human) isotope values is going to be challenging.
Future studies should take advantage of the regional distribution in δ 34 S values in the territory of Estonia, supporting this with other isotopic proxies (δ 13 C and δ 15 N from bone and dentine collagen but also δ 13 C, δ 18 O and 87 Sr/ 86 Sr from bone carbonate and tooth enamel) to help distinguish between sulfur isotopic variation caused by origin and that caused by marine influence.